---
title: "Replication code for POLITICAL SELF-CENSORSHIP IN AUTHORITARIAN STATES: THE SPATIAL-TEMPORAL DIMENSION OF TROUBLE"
output: html_document
---

```{r setup, include=FALSE}
library(RPostgreSQL)
#library(ggplot2)

library(glmmML)
library(survival)
library(bife)
library(logistf)
# library(Zelig)
# library(brglm2)
library(lme4)
# library(rpql)
library(arm)
library(stargazer)

```

## Load data

```{r cars}

rm(list=ls(all=TRUE))

# set working directory
setwd('C:/Users/cc672/Dropbox/AC/replication_files')

# load data
load("./data/anniversaries.Rda")
load("./data/focal.Rda")
load("./data/beijinguser.Rda")
```


## Regression

## Focal Main Results
```{r}

focal.dat$weekend <- ifelse(focal.dat$workday == 1, 0, 1)

# expanded political places includes: 3 focal places, party and state headquarters, and Beijing government offices
focal.dat$expandedplaces <- ifelse(focal.dat$lu501 == 1 | focal.dat$focalplaces == 1 | focal.dat$nationalgov == 1, 1, 0)

# public places includes political places, public squares, parks, etc.
focal.dat$publicplaces <- ifelse(focal.dat$lu501 == 1 | focal.dat$focalplaces == 1 | focal.dat$nationalgov == 1 | focal.dat$publicplc == 1 | focal.dat$lu402 == 1 | focal.dat$lu505 == 1, 1, 0)

# descriptive statistics
# summary(focal.dat)
# xtabs(~focalplaces + disapproving, data=focal.dat)
# xtabs(~focalplaces + political, data=focal.dat)


listVariables <- c("disapproving", "political", "focalplaces", "expandedplaces", "publicplaces", "timing")
stargazer(dplyr::select(focal.dat, listVariables), type = "text",
          covariate.labels = c("Disapproving Posts", "Political Posts", "3 Focal Places","Political Places", "Public Places", "Focal Events")
          )

# Table 2 Model 1
h1 <- glm(disapproving ~focalplaces + timing + timing*focalplaces, data=focal.dat, family = "binomial")
summary(h1)

# Table 2 Model 2
h2<- glm(political ~ focalplaces + timing + timing*focalplaces, data=focal.dat, family = "binomial")
summary(h2)

stargazer(h1, h2,
          type="text", ci.level = 0.95, star.cutoffs = c(0.05, 0.01, 0.001),
          covariate.labels = c("3 Focal Places", "Focal Times", "Focal Times*3 Focal Places","Constant"),
          dep.var.labels   =  c("Disapproving Posts", "Political Posts")
          )


# Table 2 Model 3
h1b <- glm(disapproving ~ timing + expandedplaces + timing*expandedplaces,
                data=focal.dat, family = "binomial")
summary(h1b)

# Table 2 Model 4
h2b<- glm(political ~ timing + expandedplaces + timing*expandedplaces, data=focal.dat, family = "binomial")
summary(h2b)

stargazer(h1b, h2b, type="text",
          ci.level = 0.95, star.cutoffs = c(0.05, 0.01, 0.001),
          covariate.labels = c( "Focal Times", "Expanded Political Places", "Focal Times*Expanded Political Places",
                               "Constant"),
          dep.var.labels   =  c("Disapproving Posts", "Political Posts"))


# Table 2 Model 5
h1c <- glm(disapproving ~ timing + publicplaces + timing*publicplaces, data=focal.dat, family = "binomial")
summary(h1c)


# Table 2 Model 6
h2c<- glm(political ~ timing + publicplaces + timing*publicplaces,
              data=focal.dat, family = "binomial")
summary(h2c)


stargazer(h1c, h2c,
          type="text", ci.level = 0.95, star.cutoffs = c(0.05, 0.01, 0.001),
          covariate.labels = c("Focal Times", "Public Places", "Focal Times*Public Places","Constant"),
          dep.var.labels   =  c("Disapproving Posts", "Political Posts")
          )


```

## Clustered SE
```{r}
# library(clusterSEs)
library(miceadds)

focal.dat$weekday.f <- factor(focal.dat$weekday)
# focal.dat$event.f <- factor(focal.dat$event)

# h1fe1 <- glm(disapproving ~focalplaces + timing + timing*focalplaces + weekday.f, data=focal.dat, family = "binomial")

mod1 <- miceadds::glm.cluster( data=focal.dat, formula=disapproving ~focalplaces + timing + timing*focalplaces + weekday.f,
                cluster="uid", family="binomial")

summary(mod1)

# h3fe1 <- glm(political ~focalplaces + timing + timing*focalplaces + weekday.f, data=focal.dat, family = "binomial")
# summary(h3fe1)

mod2 <- miceadds::glm.cluster( data=focal.dat, formula=political ~focalplaces + timing + timing*focalplaces + weekday.f,
                cluster="uid", family="binomial")

summary(mod2)

# h1fe2 <- glm(disapproving ~ timing + expandedplaces + timing*expandedplaces + weekday.f, data=focal.dat, family = "binomial")
# summary(h1fe2)

mod3 <- miceadds::glm.cluster( data=focal.dat, formula=disapproving ~ timing + expandedplaces + timing*expandedplaces + weekday.f,cluster="uid", family="binomial")

summary(mod3)

# h3fe2 <- glm(political ~ timing + expandedplaces + timing*expandedplaces + weekday.f, data=focal.dat, family = "binomial")
# summary(h3fe2)

mod4 <- miceadds::glm.cluster( data=focal.dat, formula=political ~ timing + expandedplaces + timing*expandedplaces + weekday.f,cluster="uid", family="binomial")

summary(mod4)

# h1fe3 <- glm(disapproving ~ timing + publicplaces + timing*publicplaces + weekday.f, data=focal.dat, family = "binomial")
# summary(h1fe3)

mod5 <- miceadds::glm.cluster( data=focal.dat, formula=disapproving ~ timing + publicplaces + timing*publicplaces + weekday.f,cluster="uid", family="binomial")

summary(mod5)

# h3fe3 <- glm(political ~ timing + publicplaces + timing*publicplaces + weekday.f, data=focal.dat, family = "binomial")
# summary(h3fe3)

mod6 <- miceadds::glm.cluster( data=focal.dat, formula=political ~ timing + publicplaces + timing*publicplaces + weekday.f,cluster="uid", family="binomial")

summary(mod6)


```



## Placebo Test 1, Residences
```{r}
# Table 3 Model 1
h1d <- glm(disapproving ~ timing + lu101 + timing*lu101, data=focal.dat, family = "binomial")
summary(h1d)


# Table 3 Model 2
h2d<- glm(political ~ timing + lu101 + timing*lu101,
              data=focal.dat, family = "binomial")
summary(h2d)

stargazer(h1d, h2d,
          type="text", ci.level = 0.95, star.cutoffs = c(0.05, 0.01, 0.001),
          covariate.labels = c("Focal Times", "Placebo Places", "Focal Times*Placebo Places","Constant"),
          dep.var.labels   =  c("Disapproving Posts", "Political Posts")
          )
```


## Placebo Test 2: anniversaries
```{r}

anniversaries.dat$weekend <- ifelse(anniversaries.dat$workday == 1, 0, 1)


anniversaries.dat$expandedplaces <- ifelse(anniversaries.dat$lu501 == 1 | anniversaries.dat$focalplaces == 1 | anniversaries.dat$nationalgov == 1, 1, 0)


anniversaries.dat$publicplaces <- ifelse(anniversaries.dat$lu501 == 1 | anniversaries.dat$focalplaces == 1 | anniversaries.dat$nationalgov == 1 | anniversaries.dat$publicplc == 1 | anniversaries.dat$lu402 == 1 | anniversaries.dat$lu505 == 1, 1, 0)

# descriptive summary
# summary(anniversaries.dat)

library(dplyr)
listVariables <- c("disapproving", "political", "focalplaces", "expandedplaces", "publicplaces", "timing")
stargazer(dplyr::select(anniversaries.dat, listVariables), type = "text",
          covariate.labels = c("Disapproving Posts", "Political Posts", "3 Focal Places","Political Places", "Public Places", "Focal Events")
          )


# Table 4 Model 1
h1 <- glm(disapproving ~focalplaces + timing + timing*focalplaces, data=anniversaries.dat, family = "binomial")
summary(h1)

# Table 4 Model 2
h2<- glm(political ~ focalplaces + timing + timing*focalplaces, data=anniversaries.dat, family = "binomial")
summary(h2)

stargazer(h1, h2,
          type="text", ci.level = 0.95, star.cutoffs = c(0.05, 0.01, 0.001),
          covariate.labels = c("3 Focal Places", "Placebo Times", "Placebo Times*3 Focal Places","Constant"),
          dep.var.labels   =  c("Disapproving Posts", "Political Posts")
          )


# Table 4 Model 3
h1b <- glm(disapproving ~ timing + expandedplaces + timing*expandedplaces,data=anniversaries.dat, family = "binomial")
summary(h1b)


# Table 4 Model 4
h2b<- glm(political ~ timing + expandedplaces + timing*expandedplaces, data=anniversaries.dat, family = "binomial")
summary(h2b)

stargazer(h1b, h2b,
          type="text", ci.level = 0.95, star.cutoffs = c(0.05, 0.01, 0.001),
          covariate.labels = c("Placebo Times", "Expanded Political Places", "Placebo Times*Expanded Political Places","Constant"),
          dep.var.labels   =  c("Disapproving Posts", "Political Posts")
          )


# Table 4 Model 5

h1c <- glm(disapproving ~ timing + publicplaces + timing*publicplaces, data=anniversaries.dat, family = "binomial")
summary(h1c)

# Table 4 Model 6

h2c<- glm(political ~ timing + publicplaces + timing*publicplaces,
              data=anniversaries.dat, family = "binomial")
summary(h2c)

stargazer(h1c, h2c,
          type="text", ci.level = 0.95, star.cutoffs = c(0.05, 0.01, 0.001),
          covariate.labels = c("Placebo Times", "Public Places",  "Placebo Times*Public Places","Constant"),
          dep.var.labels   =  c("Disapproving Posts", "Political Posts")
          )

```

## Appendix Table B.1
```{r}
# disapproving posts
model1 <- bife(critical ~ geotag | uid, model = "logit", bias_corr = "ana", data = data)
summary(model1)

# political posts
model2 <- bife(pol ~ geotag | uid, model = "logit", data = data, bias_corr = "ana")
summary(model2)

```



Note that the `echo = FALSE` parameter was added to the code chunk to prevent printing of the R code that generated the plot.
